## This code creates Table 5 and Figure 9.

# Table 5

neighbor_dist_cat <- seq(10,100,10)

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1924,]

coercive_specs <- c("conflict_coercive_NoGEO","conflict_coercive_NoFE","conflict_coercive_FULL",
                    "conflict_m_police_NoGEO","conflict_m_police_NoFE","conflict_m_police_FULL")

conflict_coercive_NoGEO <- conleyreg(conflict_count~coercive+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
conflict_coercive_NoFE <- conleyreg(conflict_count~coercive+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
conflict_coercive_FULL <- conleyreg(conflict_count~coercive+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
conflict_m_police_NoGEO <- conleyreg(conflict_count~m_police+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
conflict_m_police_NoFE <- conleyreg(conflict_count~m_police+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
conflict_m_police_FULL <- conleyreg(conflict_count~m_police+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

coercive_table <- matrix(nrow = 7,ncol=7)

for(i in 1:length(coercive_specs)){
  coercive_table[1,i+1] <- coercive_specs[i]
  model <- get(coercive_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  coercive_table[2,i+1] <- coef
  coercive_table[3,i+1] <- paste0("(",se,")")
  coercive_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  coercive_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  coercive_table[6,i+1] <- nobs(model)
  coercive_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(coercive_table)

# Figure 9 (a)

coercive_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(conflict_count~coercive_sd+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1924 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  coercive_BAGO[j,"coefficient"] <- model[2,1]
  coercive_BAGO[j,"se"] <- model[2,2]
  coercive_BAGO[j,"upper"] <- coercive_BAGO[j,"coefficient"] + (1.96*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"lower"] <- coercive_BAGO[j,"coefficient"] - (1.96*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"upper_10pc"] <- coercive_BAGO[j,"coefficient"] + (1.645*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"lower_10pc"] <- coercive_BAGO[j,"coefficient"] - (1.645*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"pch"] <- ifelse(sign(coercive_BAGO[j,"upper"])==sign(coercive_BAGO[j,"lower"]),"black",ifelse(sign(coercive_BAGO[j,"upper_10pc"])==sign(coercive_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/conflict_coercive.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,coercive_BAGO$lower/2,x.axis,coercive_BAGO$upper/2)
points(x.axis,coercive_BAGO$coefficient/2, pch=21, bg=coercive_BAGO$pch)
abline(h=0, lty=3)
dev.off()

# Figure 9 (b)

m_police_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(conflict_count~m_police_sd+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1924 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  m_police_BAGO[j,"coefficient"] <- model[2,1]
  m_police_BAGO[j,"se"] <- model[2,2]
  m_police_BAGO[j,"upper"] <- m_police_BAGO[j,"coefficient"] + (1.96*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"lower"] <- m_police_BAGO[j,"coefficient"] - (1.96*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"upper_10pc"] <- m_police_BAGO[j,"coefficient"] + (1.645*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"lower_10pc"] <- m_police_BAGO[j,"coefficient"] - (1.645*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"pch"] <- ifelse(sign(m_police_BAGO[j,"upper"])==sign(m_police_BAGO[j,"lower"]),"black",ifelse(sign(m_police_BAGO[j,"upper_10pc"])==sign(m_police_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/conflict_m_police.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,m_police_BAGO$lower/2,x.axis,m_police_BAGO$upper/2)
points(x.axis,m_police_BAGO$coefficient/2, pch=21, bg=m_police_BAGO$pch)
abline(h=0, lty=3)
dev.off()

